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^H Abstract. Physically meaningful periodic solutions to certain integrable partial dif- 

^3 ferential equations are given in terms of multi-dimensional theta functions associated 

vN to real Riemann surfaces. Typical analytical problems in the numerical evaluation of 

QJ these solutions are studied. In the case of hyperelliptic surfaces efficient algorithms exist 

QJ even for almost degenerate surfaces. This allows the numerical study of solitonic limits. 

^J For general real Riemann surfaces, the choice of a homology basis adapted to the anti- 

holomorphic involution is important for a convenient formulation of the solutions and 

^*^ smoothness conditions. Since existing algorithms for algebraic curves produce a homol- 

, I ogy basis not related to automorphisms of the curve, we study symplectic transformations 

^S^ to an adapted basis and give explicit formulae for M- curves. As examples we discuss solu- 

C^ tions of the Davey-Stewartson and the multi-component nonlinear Schrodinger equations. 
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1. Introduction 



^ The importance of Riemann surfaces for the construction of almost periodic solutions to 

K^ various integrable partial differential equations (PDEs) was realized at the beginning of the 

1970s by Novikov, Dubrovin and Its, Matveev. The latter found the Its-Matveev formula 
for the Korteweg-de Vries (KdV) equation in terms of multi-dimensional theta functions on 
CN hyperelliptic Riemann surfaces. Similar formulae were later obtained for other integrable 

f"^ PDEs as nonlinear Schrodinger (NLS) and sine-Gordon equations. For the history of the 

O topic the reader is referred to the reviews [2J and [9J. Krichever [21] showed that theta- 

functional solutions to the Kadomtsev-Petviashvili equation can be obtained on arbitrary 
Riemann surfaces. The problems of real-valuedness and smoothness of these solutions were 
solved by Dubrovin and Natanzon in [llj. 

Novikov criticized the practical relevance of theta functions since no numerical algo- 
rithms existed at the time to actually compute the found solutions. He suggested an 
effective treatment of theta functions (see, for instance, [9j) by a suitable parametrization 
of the characteristic quantities of a Riemann surface, i.e., the periods of holomorphic and 
certain meromorphic differentials on the given surface. This program is limited to genera 
smaller than 4 since so-called Schottky relations exist for higher genus between the compo- 
nents of the period matrix of a Riemann surface. The task to find such relations is known 
as the Schottky problem. This led to the famous Novikov conjecture for the Schottky 
problem that a Riemann matrix (a symmetric matrix with negative definite real part) is 
the matrix of ;B-periods of the normalized holomorphic differentials of a Riemann surface 
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if and only if Krichever's formula with this matrix yields a solution to the KP equation. 
The conjecture was finally proven by Shiota [26j. 

First plots of KP solutions appeared in [24j and via Schottky uniformizations in [4]. 
Since all compact Riemann surfaces can be defined via non-singular plane algebraic curves 
of the form 



N M 



(1.1) F(x, I/) := 2^ 2^ amnx^y^ = 0, x, ^ E C, 

n=l m=l 

with constant complex coefficients anm^ Deconinck and van Hoeij developed an approach 
to the symbolic-numerical treatment of algebraic curves. This approach is distributed as 
the algcurves package with Maple, see O [71 [8]. A purely numerical approach to real 
hyperelliptic Riemann surfaces was given in [TT^TBJ, and for general Riemann surfaces in 
[16j . For a review on computational approaches to Riemann surfaces the reader is referred 
to 0. 

In this paper we want to address typical analytical problems appearing in the numerical 
study of theta-functional solutions to integrable PDEs, and present the state of the art 
of the field by considering concrete examples. The case of hyperelliptic Riemann surfaces 
{N = 2 in ( |1.1| )) is the most accessible, since equation (1.1) can be solved explicitly for 



y, and since a basis for differentials and homology can be given a priori. Families of 
hyperelliptic curves can be conveniently parametrized by their branch points. The codes 
|14( [T5] are able to treat effectively numerically collisions of branch points, a limit in 
which certain periods of the corresponding hyperelliptic surface diverge. If the limiting 
Riemann surface has genus 0, the theta series breaks down to a finite sum which gives for 
an appropriate choice of the characteristic well known solitonic solutions to the studied 
equation. 



For solutions defined on general real algebraic curves, i.e., curves (1.1) with all anm 
real, an important point in applications are reality and smoothness conditions. These are 
conveniently formulated for a homology basis for which the ^-cycles are invariant under 
the action of the anti-holomorphic involution. However, the existing algorithms for the 
computational treatment of algebraic curves produce a basis of the homology that is in 
general not related to possible automorphisms of the curve. To implement the reality 
and smoothness requirements, a transformation to the basis for which the conditions are 
formulated has to be constructed. We study the necessary symplectic transformations and 
give explicit relations for so-called M-curves, curves with the maximal number of real ovals. 

To illustrate these concepts, we study for the first time numerically theta-functional 
solutions to integrable equations from the family of NLS equations, namely, the multi- 
component nonlinear Schrodinger equation 



dil^j d tp 



2./, / ^ 



^^■^^ iif+ a^ + MZ.^^IV'fcl K.-o, j 



, n, 



denoted by n-NLS^, where s = (si, . . . , 5^), Sk = ±1, and the (2 + l)-dimensional Davey- 
Stewartson (DS) equations, 

i i^t + i^xx -a'^^pyy + 2{^ + p l^p) ^ = 0, 
(1.3) $,, + a2$^^ + 2p|^|2^ = 0, 



where o: = i or o: = 1 and where p = ±1. Both equations (1.2) and (1.3) reduce to the 



NLS equation under certain conditions: the former obviously in the case n = 1, the latter 



if '0 is independent of the variable y and satisfies certain boundary conditions, for instance 
that $ + p |'0p tends to zero when x tends to infinity. 

Integrabihty of the NLS equation was shown by Zakharov and Shabat [3T] and algebro- 
geometric solutions were given by Its [18j. The multi-component nonlinear Schrodinger 
equation ( 1.2| ) in the case n = 2, 5 = (1,1), is called the vector NLS or Manakov system. 
Manakov [23j first examined this equation as an asymptotic model for the propagation of 
the electric field in a waveguide. Its integrability was shown for n = 2 by Zakharov and 
Schulman in [32j and for the general case in [25]. Algebro-geometric solutions to the 2-NLS^ 
equation with s = (1,1) were presented in fT2], and for the general case in [19j. The DS 
equation ( |1.3| ) was introduced in [5j to describe the evolution of a three-dimensional wave 
packet on water of finite depth. Its integrability was shown in [Ij, and solutions in terms 
of multi-dimensional theta functions on general Riemann surfaces were given in [22} [T9]. 

To ensure the correct numerical implementation of the formulae of [19], we check for 
each point in the spacetime whether certain identities for theta functions are satisfied. 
Since these identities are not used in the code, they provide a strong test for the computed 
quantities. Numerically the identities are never exactly satisfied, but to high precision. 
The code reports a warning if the residual of the test relations is larger than 10~^ which is 
well below plotting accuracy. Typically the conditions are satisfied to machine precisioijj 
In addition we compute the solutions on a numerical grid and numerically differentiate 
them. We check in this way for low genus that the solutions to n-NLS^ and DS in terms 
of multi-dimensional theta functions satisfy the respective equations to better than 10~^. 
These two completely independent tests ensure that the presented plots are showing the 
correct solutions to better than plotting accuracy. 

The paper is organized as follows: in Section 2 we recall some facts from the theory of 
multi-dimensional theta functions and the theory of real Riemann surfaces, necessary to 
give theta-functional solutions to the n-NLS^ and DS equations. In Section 3 we consider 
the hyperelliptic case and study concrete examples of low genus, also in almost degenerate 
situations. In Section 4 we consider examples of non- hyperelliptic real Riemann surfaces 
and discuss symplectic transformations needed to obtain smooth solutions. We add some 
concluding remarks in Section 5. 

2. Theta functions and real Riemann surfaces 

In this section we recall basic facts on Riemann surfaces, in particular real surfaces, and 
multi-dimensional theta functions defined on them. Solutions to the n-NLS^ and the DS 
equations in terms of theta functions will be given following [T9]. 

2.1. Theta functions. Let TZg be a compact Riemann surface of genus ^ > 0. Denote 
by (^, B) '.— (^1, . . . , Ag^ ;Bi, . . . , Bg) a canonical homology basis, and by (cji, . . . , ujg) the 
basis of holomorphic differentials normalized via 



(2.1) / ujj^2m5kj, A:,j = 1,...,^. 

JAk 

The matrix B with entries Mj^j — J^ ujj of ;B-periods of the normalized holomorphic dif- 
ferentials cjj, j = 1, . . . ,^, is symmetric and has a negative definite real part. The theta 



We work with double precision, i.e., a precision of 10 ^^; due to rounding errors this is typically 
reduced to 10"^^ to lO"^'^. 
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function with (half integer) characteristic 5 = [5i, ^2] is defined by 

(2.2) eB[5](z)= J2 exp{|(B(m + 5i),m + 5i) + (m + 5i,z + 2i^52)}, 

for any z G C^; here 5i, ^2 G {O, ^} are the vectors of the characteristic 6; (., .) denotes 
the scalar product (u, v) = J^i^i^i ^^^ ^^Y u, v G C^. The theta function 0[5](z) is even 
if the characteristic S is even, i.e., if 4 (5i, ^2) is even, and odd if the characteristic 6 is odd, 
i.e., if 4 (5i, ^2) is odd. An even characteristic is called non-singular if O[5](0) ^ 0, and an 
odd characteristic is called non-singular if the gradient VO[5](0) is non-zero. The theta 
function with characteristic is related to the theta function with zero characteristic (the 
Riemann theta function denoted by O) as follows 

(2.3) e[5](z) = e(z + 2i7T62 + Ml) exp {^(B5i, 5i) + (z + 21^52, 61)} . 

Denote by A the lattice A = {2i7rN + BM, N,M G Z^} generated by the A and B- 
periods of the normalized holomorphic differentials cjj, j = 1, . . . ,^. The complex torus 
J {Tig) — C^/A is called the Jacobian of the Riemann surface TZg. The theta function (2.2) 
has the following quasi-periodicity property with respect to the lattice A: 

e[5](z + 2i7rN + BM) 

(2.4) = e[5](z)exp{-^(BM,M) - (z, M) + 2i^((5i, N) - (52,M))} . 

For the formulation of solutions to physically relevant integrable equations in terms of 
multi-dimensional theta functions, there is typically a preferred homology basis in which 
the solution takes a simple form. Let {A^ B) and (^4, B) be arbitrary canonical homology 
basis defined on TZg^ represented here by 2^-dimensional vectors. Under the change of 
homology basis 

(c ^) (J) - ii 



where I „ n ) ^ Sp{2g,'L) is a symplectic matrix, the theta function (2.2) transforms as 



(2.6) eB[5](z) = K \/dett exp |i z* {K^)-^Bi\ %[5\{z), 

where IK = 2mA + BM and 

(2.7) B = 2i7r(2i7rC + L>B)]K-\ 

(2.8) z = (2i7r)-^]K*z, 
5i\ { A -B\ (5i\ 1^. (BA^ 



for any z G C^, where Diag(.) denotes the column vector of the diagonal entries of the 
matrix. Here k is a constant independent of z and B (the exact value of /^ is not needed 
for our purposes). 

The Abel map TZg — > JiJZg) is defined by 



fP fP 

(2.10) / := / c^, 



fpo Jpo 

for any p G TZg^ where po G TZg is the base point of the application, and where u = 
(cji, . . . , (jjg)^ is the vector of the normalized holomorphic differentials. 



Now let ka denote a local parameter near a G TZg and consider the following expansion 
of the normalized holomorphic differentials cjj, j = 1, . . . , ^, 

(2.11) LOjip) = (Vaj + Waj kaip) + • • •) dfc„(p), 

for any point p G TZg lying in a neighbourhood of a, where Vaj, Waj G C. Let us 
denote by Da (resp. D^^) the operator of the directional derivative along the vector 
Va = (Va^i, • • . ,Va^gY (rcsp. W^). According to [24j and [19], the theta function sat- 
isfies the following identities derived from Fay's identity fT3]: 

(2.12) I).Z..lne(z) = ,, + ,.^^ + ^«^®(^-^'^ 



e(z)2 

(2.13) 

^"^" e(z) +^"^" e(z) +r"^^ e(z) ^^V +2i?^ine(z) + i^2 = o, 

for any z G C^ and any distinct points a, 6 G 7^^; here the scalars qi^Ki^ i = 1,2 depend 
on the points a, b and are given by 

(2.14) qi{a,b)^DaD,\ne[S]{jl^), 

(2.15) ,.(a,6) = ^^^MI^MM 



0[^](/:: 



2 



(2^17) A'2(<1, 6) = - c; lne(/„') - Dj In (e(/„') 6(0)) - (d. Ine(Xf) - A'^o, 6) 

where 5 is a non-singular odd characteristic. 

2.2. Real Riemann surfaces. A Riemann surface TZg is called real if it admits an anti- 
holomorphic involution r : TZg -^ TZg, r^ = id. The connected components of the set of 
fixed points of the anti-involution r are called real ovals of r. We denote by 7^^(K) the set 
of fixed points. According to Harnack's inequality |17| , the number k of real ovals of a real 
Riemann surface of genus g cannot exceed ^ + 1: 0<A:<^+1. Curves with the maximal 
number A: = ^ + 1 of real ovals are called M-curves. 

The complement TZg\7lg(M.) has either one or two connected components. The curve TZg 
is called a dividing curve if TZg \7^^(M) has two components, and TZg is called non-dividing 
if TZg \ TZg(M.) is connected (notice that an M-curve is always a dividing curve). 

Example 2.1. Consider the hyperelliptic curve of genus g defined by the equation 

2g+2 

(2.18) /i2= n(^-^^)' 

where the branch points A^ G M are ordered such that Ai < . . . < A2^+2- On such a 
curve, we can define two anti-holomorphic involutions ri and T2, given respectively by 
Ti(A,/i) = (A,/i) and T2(A,/x) = (A, — /Z). Projections of real ovals of ti on the X-plane 
coincide with the intervals [A2^+2, Ai], . . . , [A2^, A2^+i]; whereas projections of real ovals of 
T2 on the X-plane coincide with the intervals [Ai, A2], . . . , [A2^+i, A2^+2]- Hence the curve 



(2.18) is an M-curve with respect to both anti-involutions ti and T2. 
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Let (v4, 0) be a basis of the homology group Hi{7lg). According to Proposition 2.2 in 
Vinnikov's paper [30j, there exists a canonical homology basis (that we call for simplicity 
'Vinnikov basis' in the following) such that 



(2.19) 



tA 
tB 







where !„ is the gxg unit matrix, and H is a block diagonal gxg matrix, defined as follows: 



1) if Uc 



+ 



/O 1 
1 



/I 



1 

1 



V 

rank(BI) = ^ + 1 — A: in both cases. 



0/ 



if TZg is dividing, 



0/ 



if TZg is non-dividing; 



2)ifK 



= 0, (i.e. the curve does not have real oval), then 

/O 1 



/O 1 
1 



V 



\ 



1 

1 0/ 



or 



1 



1 

1 



0/ 



V 

rank(H) = ^ if ^ is even, rank(BI) = ^ — 1 if ^ is odd. 

Now let us choose the canonical homology basis in Hi(JZg) satisfying (2.19), take a, 6 G 
TZg and assume that ra — a and rh — h. Denote by ^ a contour connecting the points a 
and h which does not intersect the canonical homology basis. Then the action of r on the 
generators (^, 0, t) of the relative homology group Hi(TZg^ {a, h}) is given by 



(2.20) 



where the vectors N, M G Z^ are related by (see [19j) 
(2.21) 2N + HM = 0. 




2.3. Theta-functional solutions of the n-NLS^ equation. Algebro-geometric data 



associated to smooth theta-functional solutions of the n-NLS^ equation (1.2) consist of 



{TZg, T, f,Za}^ where TZg is a compact Riemann surface of genus ^ > dividing with respect 
to an anti-holomorphic involution r, and admitting a real meromorphic function / of degree 
n + 1; here z^ ^ I^ is a non critical value of / such that the fiber f~^{za) = {ai, . . . , a^+i} 
over Za belongs to the set 7^^(R). Let us choose natural local parameters ka^ near aj given 
by the projection map /, namely, ka^ip) = f{p) — Za for any point p G TZg lying in a 
neighbourhood of aj. 

Denote by {A,B,ij) the generators of the relative homology group Hi(TZg, {a^+i,aj}). 
Let d G M^ and 6^ G M. Then the following functions '^j, j = 1, . . . , n, give smooth solutions 
of the n-NLS^ equation (fOl), see [19j, 



(2.22) 



^j{x,t) 



\Aj\e 



iO 



e(z - d + Tj 
e(z - d) 



exp {-i {Ej X - Fj t)} . 



where \Aj\ = |^2(ttn+i, ttj)p^ exp {^ (d, Mj)} . The vector Mj G Z^ is defined by the 



action of r on the relative homology group Hi (JZg, {an^i^aj}) (see (2.20)). Moreover, 
Tj = J^ cj, and the vector Z reads 



where vectors V^^^-^ and W^^^^^ are defined in (2.11). The scalars Ej^Fj are given by 
(2.23) Ej = Ki{an+i,aj), Fj = ^2(^^+1,^^) - 2 Vgi(a^+i, a^). 



k=l 



and scalars qi^Ki, z = 1,2 are defined in (2.14)-(2.17). According to [19], necessary condi- 



tions for the functions ipj in (2.22) to solve the n-NLS^ equation are the identities (2.12) 



and (2.13) with (a, 6) := (a^+i, aj). 



The signs 5i, . . . , 5^ in (1.2) are given by 

(2.24) Sj = exp {i7r(l + aj)} , 

where aj G Z denote certain intersection indices on TZg defined as follows: let CLn+i^dj G 
7^^(R) lie in a neighbourhood of a^+i and aj respectively such that /(a^+i) = fisij). 
Denote by Ij an oriented contour connecting a^+i and hj. Then 

(2.25) aj = {rlj-ij)o£j 



is the intersection index of the closed contour rij 



ij and the contour ij; this index is 



computed in the relative homology group Hi{TZg^ {^n+i, <^j})- 

In particular, it was shown in [19] that solutions of the focusing n-NLS^ equation, i.e., 
for 5 = (1, . . . , 1), are obtained when the branch points of the meromorphic function / are 
pairwise conjugate. 



2.4. Theta-functional solutions of the DS equations. Now let us introduce smooth 
solutions of the DS equations. In characteristic coordinates 



C= i^{x-\ay), T] 



1 



(x + iay), 



a 



ior 1, 
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the DS equations (1.3) take the form 



i^t + -(% + drjrj)^ + 20^ = 0, 



(2.26) 



d^drj^ 



+p^(% + a,,)l^p 



0, 



where := $ + p |'0p, p = ±1. Recall that DSl^ denotes the Davey-Stewartson equation 
when a = i (in this case ^ and 77 are both real), and DS2^ when o; = 1 (in this case ^ and 
77 are pairwise conjugate). 

In both cases, for the DSl^ and DS2^ equations, the solutions have the form p2| [T9]: 

i,e(Z-d + r) 



(2.27) 



V'(e,r?,t) = |^|e'' 



e(z - d) 



exp{-i(Gie + G'2r?-G3§)} 



(2.28) 



</.(e,r?,t) = -(lne(Z 



d))jj + -(lne(Z-d)),,+ 



Here r = J^ cj for some distinct points a, 6 G TZg^ and the vector Z is defined as 

(2.29) Z = i(KiV„e-K2V6r7 + (K?W„-KiWfe)|). 
Moreover, the scalars Gi, G2 and G3 read 

(2.30) Gi = 1^1 Ki{a, 6), G2 = /^2 Ki{h, a), 

(2.31) G3 = ^l K2{a, b) + /^l K2{b, a) + /i, 

where the scalars Ki,K2 are defined in (|2.16[) and (2.17|). As shown in [19], necessary 



(2.12) and (2.13). 



conditions for the functions tp (2.27) and (2.28) to solve the DS equations are the identities 



Algebro-geometric data associated to smooth solutions (2.27), (2.28) of the DSl^ equa- 
tion consist of {TZg^ r, a, 6, /c^, A:^}, where TZg is a compact Riemann surface of genus ^ > 0, 
dividing with respect to an anti-holomorphic involution r, a, b are two distinct points 
in 7^^(R), and ka^k^, denote local parameters near a and b respectively which satisfy 
ka{^p) — ka{p) for any p lying in a neighbourhood of a, and ki,{rp) = ki,{p) for any p 
lying in a neighbourhood of b. The remaining quantities satisfy the conditions: d G M^, 
6>,/iGM, /^2 gM\{0}, and 

(2.32) A^i = -pkj K2 q2{a, b) exp {^ (BM, M) + (r + d, M)} , 

for some ki G M, where M E Z^ is defined in ( 2.20| ). The scalar |yl| is given by 

1^1 = \iiii^2q2{ci,b)\ exp{(d, M)}, 

where the quantity q2 is defined in ( 2.15[ ). 



Algebro-geometric data associated to smooth solutions ( |2.27[ ), ( |2.28[ ) of the DS2^ equa- 
tion consist of {TZg, r, a, 6, fc^, A:^}, where TZg is a compact Riemann surface of genus ^ > 
with an anti-holomorphic involution r, a, b are two distinct points such that ra = 6, and 
/ca,/i:6 denote local parameters near a and 6 respectively which satisfy ki,(rp) = ka{p) for 
any point p lying in a neighbourhood of a. Moreover, d G iM^, 0, /i G M, /^i, /^2 ^ C \ {0} 
satisfy 7^ = /^2, and the scalar l^l is given by 

\A\^\Ki\\q2{a,b)\'/^ 

Smooth solutions of the DS2+ equation are obtained when the curve TZg is an M-curve 
with respect to r, whereas solutions to DS2~ are smooth if the associated Riemann surface 



does not have real oval with respect to r, and if there is no pseudo-real function of degree 
^ — 1 on it (i.e., function which satisfies f{rp) — —f{p)~^), see [22j. 

Remark 2.1. The symmetric structure of the DS equations \2.2^ with respect to ^ and r] 
implies that a solution i/j = ^(^,77,^) to DSl^ leads to a solution ^{—^^rj^t) of DS1~ . 

3. Hyperelliptic case 

Here we consider concrete examples for the solutions, in terms of multi-dimensional 
theta functions, to DS and n-NLS^ on hyperelliptic Riemann surfaces. We first review the 
numerical methods to visualize the solutions and discuss how the accuracy is tested. 

3.1. Computation on real hyperelliptic curves. The simplest example of algebraic 
curves are hyperelliptic curves, 

2 J rii^i {^~^i)i without branching at infinity 
1 rii^i {^~^i)i with branching at infinity 

where g is the genus of the Riemann surface, and where we have for the branch points A^ G C 
the relations A^ 7^ \j for i ^ j. If the number of finite branch points is odd, the curve is 
branched at infinity. Recall that all Riemann surfaces of genus 2 are hyperelliptic, and that 
the involution a which interchanges the sheets, a{\^ii) — (A, — //), is an automorphism on 
any hyperelliptic curve in the form ( |3.1| ). A vector of holomorphic differentials for these 
surfaces is given by (1, A, ... , A^~^)* dX/fi. For a real hyperelliptic curve, the branch points 
are either real or pairwise conjugate. As we saw in Example 2.1, if all branch points A^ are 
real and ordered such that Ai < . . . < A2^+2, the hyperelliptic curve is an M-curve with 
respect to both anti-holomorphic involutions ti and T2 defined in the example. The other 
case of interest in the context of smooth solutions to n-NLS^ and DS are real curves without 
real branch point. For the involution ti, a curve given by fi^ = Yli=i {^ ~ ^i){^ ~ ^i)i with 
A^ G C\R, i = 1, . . . , ^+1, in this case is dividing (two points whose projections onto C have 
respectively a positive and a negative imaginary part cannot be connected by a contour 
which does not cross a real oval), whereas a curve given by /i^ = — flf^i (^ ~ \){^ ~ \) 
has no real oval, and vice versa for the involution T2. 

In the following, we will only consider real hyperelliptic curves without branching at 
infinity and write the defining equation in the form /x^ = (A — ^)(A — 77) Ylf^-^{X — Ei){X — Fi). 
It is possible to introduce a convenient homology basis on the related surfaces, see Fig. [l] 
for the case 77 = C- 

The simple form of the algebraic relation between /a and A for hyperelliptic curves 
makes the generation of very efficient numerical codes possible, see, for instance, [lH [15] 
for details. These codes allow the treatment of almost degenerate Riemann surfaces, i.e., 
the case where the branch points almost collide pairwise, where the distance of the branch 
points is of the order of machine precision: \Ei — F^| ^ 10~^^. The homology basis Fig. [I] 
is adapted to this kind of degeneration. 

The Abel map J uo between two points a and h is computed in the following way: the 
sheet identified at the point a — {X{a)^fi{a)) (where we take for /x the root computed by 
Matlab) is labeled sheet 1, and at the point (A(a), — ^(a)), sheet 2. Then the ramification 
point whose projection to the A-sphere has the minimal distance to A(a) is determined. 
For simplicity we assume always that this is the point ^ in Fig. [l] (for another branch point, 
this leads to the addition of half-periods, see e.g. [2J). This means we compute J^ uj as 
j^uj — j.uj — j^ uj. The choice of a branch point as the base point of the Abel map has the 
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<^^^:-<^S?, 






s / 



Figure 1. Homology basis on real hyperelliptic curves, contours on sheet 1 
are solid, contours on sheet 2 are dashed. A-cycles are the closed contours 
entirely on sheet 1. 



advantage that a change of sheet of a point a just imphes a change of sign of the integral: 
j( {a),^{a)) ^ ^ _ j( {a)-^{a))^^ rp^ computc the integral J^(jO, one has to analytically 

continue ji on the connecting line between A(a) and ^ onto the A-sphere. Whereas the root 
/i is not supposed to have any branching on the considered path, the square root in Matlab 
is branched on the negative real axis. To analytically continue fi on the path [A(a),^], we 
compute the Matlab root at some \j G [A(a),^], j = 0, . . . , A/'c and analytically continue 
it starting from ii{a) by demanding that |/i(Aj+i) — /i(Aj)| < |/i(Aj+i) + /i(Aj)|. The so 
defined sheets will be denoted here and in the following by numbers, i.e., a point on sheet 
1 with projection A(a) into the base is denoted by {\{a)Y^\ 

Thus the computation of the Abel map is reduced to the computation of line integrals 
on the connecting line between A (a) and ^ in the complex A-plane. For the numerical 
computation of such integrals we use Clenshaw-Curtis integration (see, for instance, [27]): 
to compute an integral J-,^ h{x) dx, this algorithm samples the integrand on the A/'c + 1 
Chebyshev collocation points Xj — cos(J7r/A'c), j — 0, ...,A'c. The integral is approxi- 
mated as the sum: j_-yh{x)dx ^ J2j=o^j ^^(^j) (^^^ l22j on how to obtain the weights 
Wj). It can be shown that the convergence of the integral is exponential for analytic func- 
tions h as the ones considered here. To compute the Abel map, one uses the transformation 
A -^ A(a)(l + x)/2 + ^(1 — x)/2, to the Clenshaw-Curtis integration variable. The same 
procedure is then carried out for the integral from ^ to b. 

The theta functions are approximated as in [Mj as a sum. 



(3.2) eB[5](z) 



Ne 

E 

'mi=-Ne 



Ne 



exp {^(B(m + 5i), m + 5i) + (m + 5i, z + 21^62)} . 



-Ne 



The periodicity properties of the theta function (2.4) make it possible to write z = zq + 
2i7rN + BM for some N,M G Z^, where zq = 2i7ra + B/3 with ai,/3i G ] - j ^] for 
i = 1, . . . ,^. The value of A"^ is determined by the condition that all terms in (2.2) with 
Itti^I > No are smaller than machine precision, which is controlled by the largest eigenvalue 
of the real part of the Riemann matrix (the one with minimal absolute value since the real 
part is negative definite), see [T4t[T6i|. 



11 



To control the accuracy of the numerical solutions, we use essentially two approaches. 



First we check the theta identity (2.13), which is the underlying reason for the studied 
functions being solutions to n-NLS^ and DS, at each point in the spacetime. This test 
requires the computation of theta derivatives not needed in the solution (which slightly 
reduces the efficiency of the code since additional quantities are computed), but provides 



an immediate check whether the solution satisfies (2.13) with the required accuracy. Since 
this identity is not implemented in the code, it provides a strong test. This ensures that all 
quantities entering the solution are computed with the necessary precision. In addition, the 
solutions are computed on Chebyshev collocation points (see, for instance, [27J) for each 
of the physical variables. This can be used for an expansion of the computed solution in 
terms of Chebyshev polynomials, a so-called spectral method having in practice exponential 
convergence for analytic functions as the ones considered here. Since the derivatives of the 
Chebyshev polynomials can be expressed linearly in terms of Chebyshev polynomials, a 
derivative acts on the space of polynomials via a so called differentiation matrix. With these 
standard Chebyshev differentiation matrices (see [27J), the solution can be numerically 
differentiated. The computed derivatives allow to check with which numerical precision 
the PDE is satisfied by a numerical solution. With these two independent tests, we ensure 
that the shown solutions are correct to much better than plotting accuracy (the code 
reports a warning if the above tests are not satisfied to better than 10~^). 



3.2. Solutions to the DS equations. The elliptic solutions are the well known travelling 
wave solutions and will not be discussed here. The simplest examples we will consider for 
the DS solutions are given on hyperelliptic curves of genus 2. As we saw in Section 2.4, 
for DSl^ reality and smoothness conditions imply that the branch points of the curve are 
either all real (M-curve) or all pairwise conjugate (dividing curve). The points a and b 
must project to real points on the A-sphere and must be stable under the anti-holomorphic 
involution r (we use here r = ti, as defined in Example 2.1, except for DS2~). For 
DS2^, we have ra — h where the projection of a onto the A-sphere is the conjugate of the 
projection of h. For DS2^ the curve must have only real branch points (M-curve), whereas 
for DS2~ it must have no real oval. 

For DS we will mainly give plots for fixed time since for low genus, the solution is 
essentially travelling in one direction. For higher genus, we show a more interesting time 
dependence in Fig. [9| 

We first consider the defocusing variants, DSl^ and DS2+ on M-curves. In genus 2 
we study the family of curves with the branch points —2, —1, 0, e, 2, 2 + e for e = 1 and 
6 = 10~^^. In the former case the solutions will be periodic in the (x, ?/)-plane, in the latter 
almost solitonic since the Riemann surface is almost degenerate (in the limit e ^ the 
surface degenerates to a surface of genus 0; the resulting solutions are discussed in more 
detail in [20j). To obtain non-trivial solutions in the solitonic limit, we use d = 2 [o ol ^^ 
all examples. In Fig. [2] it can be seen that these are in fact dark solitons^ i.e., the solutions 
tend asymptotically to a non-zero constant and the solitons thus represent 'shadows' on a 
background of light. The well known features from soliton collisions for (l+l)-dimensional 
integrable equations, namely, the propagation without change of shape, and the unchanged 
shape and phase shift after the collision, can be seen here in the (x, y)-plane. 

The corresponding solutions to DS2+ can be seen in Fig. [3| We only show the square 
modulus of the solution here for simplicity. For the real and the imaginary part of such a 
solution for the DSl~-case, see Fig. ^ Because of remark 2.1 all solutions shown for DSl^ 
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on M-curves are after the change of coordinate ^ ^ — C solutions to DSl" 
DS1~ solutions on M-curves will not be presented here. 



For this reason 





Figure 2. Solution (2.27) to the DSl^ equation att — {) on a hyperelliptic 
curve of genus 2 with branch points —2, —1, 0, e, 2, 2 + e and a = (— 1.9)^"'^^; 
b = ( — 1.1)^^) for e — 1 on the left and e = lO"-*^^; the almost solitonic limit, 
on the right. 





Figure 3. Solution {2.2'^ to the DS2^ equation att — Q on a hyperelliptic 
curve of genus 2 with branch points —2, — 1, 0, e, 2, 2 + e and a — (— 1.5 + 
2i)(i), b = (-1.5 - 2i)(2) for e^l on the left and e = 10"^^ the almost 
solitonic limit, on the right. 

In the same way one can study, on a genus 4 hyperelliptic curve, the formation of 
the dark 4-soliton for these two equations. We consider the curve with branch points 
-4, -3,-2, -2 + 6,0,6,2, 2 + 6,4,4 + e for 6 = 1 and e = 10"^°, and use d = | [qooo]^- 
The DSl^ solutions for this curve can be seen in Fig. |4J The corresponding solutions to 
DS2+ is shown in Fig. [5| 

Solutions to the focusing variants of these equations can be obtained on hyperelliptic 
curves with pairwise conjugate branch points. For such solutions the solitonic limit cannot 
be obtained as above since the quotient of theta functions in ( |2.27[ ) tends to a constant 
in this case. To obtain the well-known bright solitons (solutions tend to zero at spatial 
infinity) in this way, the hyperelliptic curve has to be completely degenerated (all branch 
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Figure 4. Solution (2.27) to the DSl^ equation att — {) on a hyperelliptic 
curve of genus 4 ^ith branch points —4, —3, —2, —2 + e, 0, e, 2, 2 + e, 4, 4 + 6 
and a = {-3.9)^^\ b = (-S.l)^^) for e = 1 on the left and e = IQ-^^ the 
almost solitonic limit, on the right. 





Figure 5. Solution (2.27) to the DS2^ equation att — Q on a hyperelliptic 
curve of genus 4 ^ith branch points —4, —3, —2, —2 + e, 0, e, 2, 2 + e, 4, 4 + 6 
anda= (-1.5 + 2i)(i), 6= (-1.5-2i)(i) fore^l on the left and e = 10" 
the almost solitonic limit, on the right. 
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points must collide pairwise to double points) which leads to limits of the form '0/0' in 



the expression for the solution (2.27) which are not convenient for a numerical treatment; 



see ^20j for an analytic discussion. Therefore we only consider non-degenerate hyperelliptic 
curves here. To obtain smooth solutions, we use d = 0. A solution in genus 2 of the DS1~ 
equation is studied on the curve with the branch points — 2 zb i, — 1 ± i, 1 ± i in Fig. |6| 

A typical example of a DS1~ solution on a hyperelliptic curve of genus 4 with branch 
points — 2 zh i, — 1 ± i, zbi, 1 zh i, 2 zb i is shown in Fig. [7[ 

Smooth solutions to DS2~ can be obtained on Riemann surfaces without real oval for 
points a and b satisfying ra = b. As mentioned above, hyperelliptic curves of the form 
fi^ — — Y\iLi (^ ~ \) with pairwise conjugate branch points have no real oval for the 
standard involution ri as defined in Example 2.1. On the other hand, surfaces defined 
by the algebraic equation fj? — Y\iLi {^ ~ \) have no real oval for the involution T2 (see 
Example 2.1). We will consider here T2 for the same curves as for DS1~. An example for 
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Figure 6. Solution (2.27) to the DS1~ equation att — {) on a hyperelliptic 
curve of genus 2 with branch points — 2 ± i^ — 1 ± i^ 1 zb i and a = (— 4)(-^)^ 
b = (— 3)(^). The square modulus of the solution is shown on the left, real 
and imaginary parts on the right. 




Figure 7. Solution 12.21 ) to the DSl equation att — on a hyperelliptic 
curve of genus 4 ^ith branch points — 2 zb i^ — 1 zb i^ zbi^ 1 i i; 2 ± i and 
a = (-4)(i), 6=(-3)(2). 



genus 2 can be seen in Fig. [sj An example for a DS2 
Fig-i 



solution of genus 4 can be seen in 



3.3. Solutions to the n-NLS^ equation. A straightforward way to obtain solutions 



(2.22) to the n-NLS^ equation is given on an (n + l)-sheeted branched covering of the 
complex plane, an approach that will be studied in more detail in the next section. As can 
be seen from the proof of Theorem 4.1 in [19j, the crucial point in the construction of these 
solutions is the fact that J2k=i ^^fc ~ ^- This implies that it is also possible to construct 
theta-functional n-NLS^ solutions on hyperelliptic surfaces by introducing constants ^k via 
J2k=i ^k^ak = in the following corollary of Theorem 4.1 in [19j: 
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Figure 8. Solution to the DS2 equation att — {) on a hyperelliptic curve 
of genus 2 with branch points — 2 zb i^ — 1 zb i^ 1 ± i and a — ( — 1.5 + 2i)(-'^)^ 
6=(-1.5-2i)(2). 



t=o 



t=5 




t=10 



t=20 




Figure 9. Solution to the DS2 equation for several values of t on a hy- 
perelliptic curve of genus 4 ^ith branch points — 2zbi^ — lii; ^h ^^h 2±i 
and a = (-1.5 + 2i)(i), b = (-1.5 - 2i)(2). 
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Corollary 3.1. Let TZg he a real hyperelliptic curve of genus g > and denote by r an 



anti-holomorphic involution. Choose the canonical homology basis which satisfies (2.19). 
Take n > g and let ai, . . . , a^+i G TlgO^) not ramification points having distinct projection 
X{aj), j = 1, . . . , n + 1^ onto the X-sphere. Denote by ij an oriented contour between a^+i 
and aj which does not intersect cycles of the canonical homology basis. Let dR G M^; 



T G Z9, and define d = di^ + f{diag{m) - 2T). Take 6 el 
be arbitr ary co nstants with 7n+i = 1- Put s = (sign{'yi) si^ . 



and let 7^+1,..., 7^ G 



where Sj is 



l,...,g, are defined by YJk=i Ik^au = 0- ^^^n 



given in (2.24), ^nd the scalars 7j, j 

1, . . . , ri; give solutions of the n-NLS^ equation (1.2) 

e(z 



the following functions ipj , j 
(3.3) i,jix,t)^\j,\^/^\Aj 



Jd 



d + r, 



where \Aj\ — \q2{an+i,aj 



e(z-d 

1/2 exp { I (d, M^) } . Here Z = i V« 



exp{-i {Ej X — Fjt)}, 



vectors V^^^-^ and W^^^-^ were introduced in (2.11), and Vj = /^,^- The scalars Ej^Fj 



are given by 



X + i^an+i ^; where the 



Ej = Ki{an+i,aj), Fj = ^2(^^+1,^^) - 2^7^ gi(a^+i, a/e), 



k=i 



where qi^ Kj for z = 1,2 are defined in (2.14)-(2.17). If TZg is dividing and if d ^ 



functions (3.3) give smooth solutions of n-NLS^ . 



As an example we consider, as for DS in genus 2, the family of curves with the branch 
points —2, —1, 0, 6, 2, 2 + e for e = 1 and e = 10~^^. In the former case the solutions will be 
periodic in the (x,t)-plane, in the latter almost solitonic. To obtain non-trivial solutions 
in the solitonic limit, we use d = ^ [q q] in all examples. 

In Fig. [10] we show the case ai = (-1.9)(^), a2 = (-1.1)(^) and as = (-1.8)(^), which 
leads to a solution of 2-NLS^ with s — (—1,-1). Interchanging a2 and a3 in the above 
example, we obtain a solution to 2-NLS^ with 5 = (1, —1) in Fig. 
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. A ' /, 



^^^/////////^^. 




Figure 10. Solution p.jy to the 2-NLS^ equation with s = (-1,-1) on 
a hyperelliptic curve of genus 2 with branch points —2, —1, 0, e, 2, 2 + e and 
ai = (-i.9)(i)^ a2 = (-1.1)(^) and as = (-1.8)(^) /or e = 1 on the left and 
6=10"-*^^; the almost solitonic limit, on the right. 

Solutions of 4-NLS^ can be studied in the same way on the hyperelliptic curve of genus 4 
with branch points —4, —3, —2, —2 + e, 0, e, 2, 2 + e, 4, 4 + e. We use d = H ] ^^^ ^^^ 
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points ai = (-3.9)(i), as = (-3.7)(i), as = (-3.5)(i), a4 = (-3.3)(^) and as = (-3.1)^. 
Since the vectors V^^. and W^^. are very similar in this case, the same is true for the 
functions '0j. Therefore, we wih only show the square modulus of the first component '0i 
in Fig. 12 for s — (1, —1, 1, —1) on the left. Interchanging a4 and a^ in this case, one gets 
a solution to 4-NLS^ with s — (—1, 1, —1, —1) which can be seen on the right of Fig. 12 











Figure 11. Solution ( fg.gp to the 2-NLS^ equation with s = (1,-1) on a 
hyperelliptic curve of genus 2 with branch points —2, — 1,0, e, 2,2 + e and 
ai = (-1.9)(i), as = (-1.8)(i) and as = (-l.l)^^) for e = 1 on the left and 
6=10"-*^^; the almost solitonic limit, on the right. 





Figure 12. Solution to the 4-NLS^ equation on a hyperelliptic curve of 
genus 4 ^^th branch points —4, —3, —2, —2 + e, 0, e, 2, 2 + e, 4, 4 + e and 
e = 1 /or ai = (-3.9)W, as = (-3.7)^, as = (-3.5)(i\ a4 = (-3.3)(i) 
and as = (— 3.1)^-'^^ which leads to s — (1,-1,1,-1)^ on the left, and 
for ai = (-3.9)W, as = (-3.7)^, as = (-3.5)(i), a4 = (-3.1)(i) and 
as = (— 3.3)^-'^^ which leads to s — (—1, 1, —1, —1) on the right. 



The almost solitonic limit e = 10~^^ produces well known solitonic patterns as shown for 
instance for the DS equation in the previous subsection. 

Hyperelliptic solutions to the n-NLS^ equation with all signs Sj satisfying 5j = 1, can 
be constructed on a curve without real branch points. To obtain smooth solutions, we use 
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d = 0. A solution of the 2-NLS^ equation is studied on the curve of genus 2 with the 



branch points —2 zt i, — 1 zb i, 1 =b i in Fig. 13 




^i£„ Ml\Will/^£, 



-///£ . 



Figure 13. Solution to the 2-NLS^ equation with s = (1, 1) on a hyperellip- 
tic curve of genus 2 with branch points — 2zbi^ — lii; lii andai = (—1.9)^-'^^ 
a2 = (-1.8)(2) and as = (-1.1)^^^. 



A typical example for a hyperelliptic 4-NLS^ solution with 5 = (1, 1, 1, 1) can be obtained 



on a curve of genus 4 with branch points — 2 ± i, — 1 zb i, ±i, 1 zb i, 2 zb i, as shown in Fig. 14 





Figure 14. Solution to the 4-NLS^ equation with s = (1, 1, 1, 1) on a hy- 
perelliptic curve of genus 4 ^ith branch points — 2zbi^ — lii; ii; lii; 2zbi 

and ai = (-3.9)^^\ as = (-3.7)(2)^ as = (-3.5)(i\ a^ = (-3.3)(2) and 
a5 = (-3.1)«. 
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4. General real algebraic curves 

The quantities entering theta-functional solutions of the DS and n-NLS^ equations are 
related to compact Riemann surfaces. Since all compact Riemann surfaces can be defined 
via compactified non-singular algebraic curves, convenient computational approaches as 
|6l E] and [16] are based on algebraic curves: differentials, homology basis and periods of 
the Riemann surface can be obtained in an algorithmic way. We refer the reader to the 
cited literature for details. The identification of the sheets of the covering defined by the 



algebraic curve (1.1) via the projection map (x^y) \-^ x, is done, as in the hyperelliptic 
case, by analytic continuation of the roots y^, i = 1, . . . , A^ for some non-critical point x^ 
on the x-sphere, along a set of contours specified in [16]. In the context of real algebraic 
curves for which solutions of n-NLS^ and DS are discussed here, an additional problem is 
to establish the action of the anti-holomorphic involution r on points on different sheets. 
A typical problem is to find points a G TZg and b G TZg with the same projection onto the 
x-sphere such that ra = 6; here r is defined via ra = (x(a),y(a)). To this end, the roots 
y^, z = 1, . . . , A^, identified at x = X5, are analytically continued to the points projecting 
to x(a) on the x-sphere. It is then established which pairs of points in the different sheets 
satisfy ra = b. 

In contrast to the hyperelliptic curves of the previous section, it is not possible for 
general curves to introduce a priori a basis of the homology. Thus the cited codes use 
an algorithm by Tretkoff and Tretkoff [28j which produces a homology basis for a given 
branching structure of the covering which is in general not adapted to possible automor- 
phisms of the curve. In the context of theta-functional solutions to integrable PDEs one 
is often interested in real curves. As discussed in [19], the Vinnikov basis (i.e., the canon- 



ical homology basis which satisfies (2.19[)) is convenient in this context. Since solutions 



and smoothness conditions for n-NLS^ and DS equations are formulated in this basis, a 
symplectic transformation relating the computed basis to the Vinnikov basis needs to be 
worked out. This transformation is discussed in the present section and will be applied to 
examples of real algebraic curves. 

4.1. Symplectic transformation. Let TZg be a real compact Riemann surface of genus 
g and r an anti-holomorphic involution defined on it. Let (z/i, . . . , Ug) be a basis of holo- 
morphic differentials such that 



(4.1) r'^Uj^Uj, j = l,...,^, 

where r* is the action of r lifted to the space of holomorphic differentials: T*a;(p) = (^(rp) 
for any p G TZg. For an arbitrary canonical homology basis (A,B)^ let us denote by P4 
and Ps the matrices of A and ;B-periods of the differentials Uj: 

(4.2) {PA)ij = / ^j, {PB)ij^ / i^j, i,j = l,...,^. 

JA^ JB^ 



In what follows {A^B) denotes the Vinnikov basis. From (4.1) and (2.19) we deduce the 



action of the complex conjugation on the matrices P4 and Ps- 

(4.3) {PA)^J e M, 

(4.4) Pb--Pb^ HP4. 

Denote by {A^ B) the homology basis on TZg produced by the Tretkoff- Tretkoff algo- 



rithm. From the symplectic transformation (2.5) we obtain the following transformation 
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law between the matrices P^, P^ and P^, Ps defined in (4.2): 



(4.5) 



A B 



C D \P, 



A 



B 



Pa 
Pb 



Therefore, by (4.3) one gets 

(4.6) A Re (P4) + 5 Re (Pg) = P4 

(4.7) A Im (P4) +Blm {Pg) = 0, 



and by (4.4) 



(4.8) 
(4.9) 



CRe{P^)+DRe{P^)^-mPA 
CIm (P^) + Dim (P^) = Im (Pg) . 



According to (4.6), the matrix A Re (P2)+P Re (Pg) is invertible, since the matrix P_4 of 
^-periods of a basis of holomorphic differentials is always invertible (see, for instance, [3j). 
Moreover, it is well known that the Riemann matrix B = 2'n: Pb {Pa)~^ has a (negative) 
definite re al p art, which is equal to — 27rIm(Pe) Im((P_4)~-^) for the real matrix P_4 here. 
Then, by (4.9) the matrix CIm {Pj) + Dim (Pg) is also invertible. 



Lemma 4.1. The matrices A,B,C,D G Mg{Z) solving U^-{4-9) satisfy: 
A^ = Im {Pg) [C Im {P^) +DIm {Pg)]~^ 



(4.10) 
(4.11) 

(4.12) 
(4.13) 



B* = -Im (P_|) [C Im {Pj} +DIm (Pg)] ^ 

& = ]- A*H - Re (Pg) [A Re {Pj) + B Re (Pg)] "^ 

D* = i P*H + Re {P^) [A Re (P^) + B Re (Pg)] "^ 



(A B 
Proof. Recall that symplectic matrices M — [ „ n ) ^ Sp{2g, Z) are characterized by 

(4.14) A^D-C^B = h 

(4.15) 

(4.16) 



'^g, 



A^C = &A, 
D^B = B^D. 



Multiplying equality (4.7) from the left by the matrix C*, we deduce from (4.14) and (4.15) 
that: 

C^A Im (P^) + C*P Im (Pg) = 

C^A Im (P_^) + {A^D - Ig) Im (Pg) = 

A' CIm {P^) + A' Dim (Pg) = Im (Pg) , 



which lea ds to (4.10). Equality ( |4.11[) c an be checked analogously with (4.14) a nd (4.16 ). 
To prove (4.12), multiply equality (4.8) from the left by the matrix A*. Using (4.14) and 
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(4.15) one gets: 



A*CRe (P^) + A*DRe {Pg) = -A^MPa 

C^ARe (P^) + {Ig + C^B) Re (Pg) = ^^*HP^ 

C* (A Re (P_4) + PRe (Pg)) = ^^*HP^ - Re (Pg) , 
which by (4.6) leads to ( |4.12[ ). Identity ( |4.13[ ) can be proved analogously. 



n 



Remark 4.1. Lemma 4.1 implies that it is sufficient to know the matrices A and B (or 
C and D) to determine the symplectic matrix in (4.5). In practice, this means that a 



convenient ansatz for one of the matrices has to be found. The others then follow from the 
relations in Lemma |4?T1 



Thus to construct these matrices one ffist checks which of the matrices Re {P^ , Re {P^ , 
Im {P^-> Im (P^) are invertible. This way a matrix can be identified (e.g. A) in terms of 
which the others can be expressed. The task is thus reduced to provide an ansatz for this 
matrix such that the others will have entire components. We illustrate this approach at 
the example of the Trott curve below. 

Proposition 4.1. Let {A^B) he the canonical homology basis obtained with the Tretkoff- 
Tretkoff algorithm; we denote with a tilde the quantities expressed in this basis. Under the 



change of homology basis (S.5i); solutions ofn-NLS^ and DS equations given in (2.22) and 



(2.21), respectively, which are expressed in the basis satisfying (2.19), transform as follows: 
the vector d appearing in the solutions becomes (2i7r)~-'^ K^ d where K = 2mA -]- BM, and 
the theta function O = Ob with zero characteristic, transforms to the theta function O^ [S\ 
with characteristic 5 — [5i,52] given by 



(4.17) 

(4.18) 
where 
(4.19) 



~5i = j diag (p*HP - 2 Re {Pj) M-^Im (p_^ 
S2 = \ diag (a^MA -2 Re (P^) M'^/m (p*^ 



Im[pi)Re{P_^)-Im(p'^)Re{P^). 



Moreover, the real constant h appearing in \2.28 ) and (2.3 1\) becomes h + h where Im{h) 
is given by 



(4.20) 



Im(h) = - In 



%[m+r) 



Im{G3 



with Z = i(Wa - W5), and the vectors N,M defined in 12. 2C^ become A^N + C^M and 
B^N + D^Nl respectively. 

Proof. Under the change of the canonical homology basis (2.5), the vector uj = (cji, . . . ^UgY 
of normalized holomorphic differentials transforms as 

(4.21) uj = 2i7T{K^)-^ 



u, 



where K = 2mA + 5B. According to the transformation law (2.6) of theta functions 



it can be checked after straightforward calculations, that under this change of homology 
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basis, quantities (2.15)-( [2.17| ) transform as: 

(4.22) q2{a,b) = q2{a,b) exp {-f* (t*)-iBf } , 

(4.23) 



1 - 



qi{a,b) ^ qi{a,b) + ^ViiK'r'BVt, 



(4.24) 



Ki{a,b)^Ki{a,b) + ^ ( V^ (K*)-^Bf + r* (K*)-^B V, 



(4.25) K2{a,b) = k2{a,b)-^ (w^^*)"'^? + r* (K*)-!^ W„) - V* (K*)-!^ V„. 



We deduce that solutions of the n-NLS* and DS equations given in (2.22) and ( |2.27 ), re- 
spectively, transform as follows: the vector d becomes (2i7r)~^ W d, and the theta function 
© = 6b with zero characteristic, transforms to the theta function 0b[5] with character- 
istic 6. To compute the vectors of the characteristic 6 we consider the inversion of the 



symplectic matrix in (2.5) which leads to 






Since the characteristic used in [19j to construct solutions (2.22), (2.27) of n-NLS* and DS 
is zero, we get with (2.9) 



Diag ( ^i^ 



(note that D^B and C^A are symmetric matrices). Substituting (4.10) and (4.11) in (4.12) 
(resp. ( |4.13| )), it can be checked that 



with 



C^A = - (A^MA - 2 Re (P^) M"^ Im f P| 
D*B = - (b*MB - 2 Re (P^) M"^ Im fp^) ) , 

M = Im [P^) Re (P^) - Im (p}) Re (P^) . 



Moreover, the real constant h appearing in the solutions ( |2.27P -(2.28) of the Davey- 
Stewartson equations becomes h + h, where h is given by 



(4.26) 



h = -\i {KT'BVa - yl o^T B\b- 



Notice that the construction of the so lutions (2.27) given in [19j allows to express the 
imaginary part of the constant h (4.26) in terms of the characteristic S. Namely, since the 
reality condition 

(4.27) V* = P^ 

is satisfied for the Vinnikov basis, where the function ip*{^,ri,t) reads 

q2{a,b) e(Z-d-r 



(4.28) 



4^*{^,V,t) = -K1K2 



A e(Z-d) 



exp{i(Gie + G'2r?-G3|)}, 
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it also holds for the computed basis. Therefore, putting ^ 

3ssed in the c 

e^[^](z + f) 



V 



0, t = 2, d = 0, and 



taking the modulus of each term in ( |4.27[ ) expressed in the computed basis, one gets: 



■exp{-i(G3 + /i)} 



where Z = i(Wa - W5). We deduce that 



emm 



exp{i{G3 + h)} 



(4.29) 



lm{h) 



1 



In 



e.[5]iZ + r) 



%m^ 



-MGs). 



D 

= 0, the vectors of 

characteristic (^.-Z?j and (4-18) do not depend explicitly on the symplectic matrix appearing 
in the change of homology basis and are uniquely defined by: 



Remark 4.2. In the case where the spectral curve is an M-curve, i.e. 



(4.30) Si ^ - dtag I Re {P^) 

(4.31) ^2 = ^ dzag Ue (P^) 



Im(pi)Re{P_^)-Im(pURe{P^) 



Im 



p'a 



Im 



Pl)Re{Pj)-Im(pOiRe{P^) 



Im 



(-1) 



It would be possible to compute the theta-functional solutions in the Vinnikov basis 
once the symplectic transformation between this basis and the basis determined by the 
code is known. However, since this symplectic transformation is not unique, the found 
Vinnikov basis leads in general to a Riemann matrix for which the theta series converges 



only slowly, i.e., the value Nq in (3.2) has to be chosen very large. To avoid this problem, 
we compute the theta function always in the typically more convenient Tretkoff-Tretkoff 



basis with the characteristic of the theta functions given by (4. 17)- (4. 19). 



4.2. Trott curve. The Trott curve [29j given by the algebraic equation 



(4.32) 



144 (x^ + /) - 225 (x^ + r ) + 350 x^y^ + 81 = 



is an M-curve with respect to the anti-holomorphic involution r defined by t(x, y) = (x^y), 
and is of genus 3. Moreover, this curve has real branch points only (and 28 real bitangents, 
namely, tangents to the curve in two places). Our computed matrices of A and ;B-periods 
reaqj 

0.0235i 0.0138i 0.0138i 




-0.0315 



0.0277i 





0.0250 



-0.0315 + 0.02351 

-0.0235i 



0.0138i 

-0.025 + 0.0277i 

0.0138i 



-0.0250 + 0.01381 
0.0250 
0.0138i 



The Trott curve being an M-curve, the vectors of the characteristic 6 satisfy (4.30) and 
([liTt, which leads to 5 = ^[??g]^ 



For the ease of representation we only give 4 digits here, though at least 12 digits are known for these 



quantities. 



24 



C. KALLA AND C. KLEIN 



A possible choice of a symplectic transformation bringing the computed basis to the 
Vinnikov basis is: 



A = 




B = 



1 








-1 









c = 




D = 




Note that the matrices A^B^C^D are not unique since the action (2.19) of the anti- 



holomorphic involution on the basic cycles allows for permutations of ^j-cycles for in- 
stance. These matrices can be computed as follows. Since the Trott curve is an M-curve, 



one has H = 0. Moreover, the matrix Im(P^) being invertible here, by (4.7) one gets: 

(4.33) B^-Alm{P^){lm{P^))-\ 
With ( |4.14[ ) and ( |4.15| ) it follows that 

(4.34) A' [D + Clm (P^) (im (Pg))"') = Is- 

The computed matrix Im {P^ (im (P^)) being (within numerical precision) equal to 



Im(P^)(lm(P^))-^= I 






and with C, P> E 7W3(Z), we get from ( |4.34D that deiA = 1. Since A G M'jj T), th e 
condition deiA — 1 implies A G G/3(Z). For any A G G/3(Z), one can see from ( |4.33 ), 
(|4.12|) and (|4.13|) that B,C,D G M^{Z), and therefore that the matrices A,B,C,D give 



a solution of (4.6)-(4.9). The choice A = I3 leads to the above matrices. 



The Trott curve has real fibers and can thus be used to construct solutions to the 3-NLS 
equation via the projection map / : (x, y) 1-^ x, which is a real meromorphic function of 
degree 4 on the curve. We consider the points on the curve stable with respect to r and 
projecting to the point with x = 0.1 in the x-sphere, and choose d = 0. The corresponding 



solution to the 3-NLS equation can be seen in Fig. 15 



A solution to the DSl^ equation on this curve can be constructed for points a and h 
stable with respect to the involution r. The solution for a — {—{).2)^^\ b = (0.2) (^^ and 
the choice d = can be seen in Fig. 16 Note that in accordance with Remark 2.1, one 
would obtain a solution of DS1~ for the choice a = (— 0.2)(-'^) and b = (0.2)(-'^). 

Similarly, a solution to the DS2+ equation can be obtained for points a and b subject 
to ra = b. For a = (0.1 + i)(^) and b = (0.1 - i)(^) we get Fig. Kjl 



4.3. Dividing curves without real branch point. We consider 
the curve given by the equation 



(4.35) 



30x^ - 61x^^ + 41y^x'^ 



43x^ - lly^x + 42xy + / - 11^^ + 9 = 



which was studied in [lOj and [30]. It is a genus 3 curve, dividing with respect to the anti- 
holomorphic involution r, without real branch point. This curve admits two real ovals. In 
this case the matrix H has the form 
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Figure 15. Solution {2.2^ to the 3-NLS^ equation on the Trott curve 
for the points with x = 0.1 on the x-sphere. The sheets are identified at 
the points projecting to x ^ -1.0129, (0.9582i, -0.9582i, 0.1146i, -0.1146i). 
The vector of signs equals 5 = (1, — 1, — 1) from top to bottom. 




Figure 16. Solution to the DSl^ equation on the Trott curve for the points 
a = (-0.2)(i) and b = (0.2)(2) att^O. 

The period matrices computed by the code read 



0.2721 - 0.0977i 


-0.3193 + 0.1914i 


-1.0668 + 0.4293i 


0.2721 + 0.0977i 


-0.3193 - 0.3341i 


-1.0668 -0.4316i 


0.2721 - 0.0977i 


0.4676 - 0.3341i 


0.7992 - 0.4316i 
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Figure 17. Solution to the DS2^ equation on the Trott curve for the points 
a = (0.1 + i)(i) and b = (0.1 - i)^^) att^O. 



0.2721 - 0.2932i 


-0.3193 + 0.33411 


-1.0668 + 0.43161 


0.2721 + 0.29321 


-0.3193 - 0.7169i 


-1.0668-1.29031 


0.2721 - 0.0977i 


0.4676 + 0.19141 


0.7992 + 0.42931 



Aft er s o me c alculations, one finds that the following matrices A, B, C, D provide a solution 

of ra-ra: 



1 


2 


-1 


2 


-1 








2 


-1 



.B 




C 



From (4.17) and (4.18) one gets for the characteristic: 6 




The curve (4.35 )nas real fibers and can thus be used to construct solutions to the 
focusing 3-NLS equation. We consider the points on the curve with x = 2.5 and stable 
with respect to r, and we choose d = 0. The corresponding solution to the focusing 3-NLS 



equation can be seen In Fig. [18 

A solution to the DS1~ equation can be constructed by choosing the points a = (— 4)^-'^) 
and 6= (-S)^^) see Fig. 
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4.4. Fermat curve. The Fermat curves 

(4.36) y^ + x"" + 1 = 0, n > 2, n even, 

are real curves without real oval with respect to r. We consider here the curve with n — A 
that has genus 3. The matrix H has the form 




27 






Figure 18. Solution to the 3-NLS^ equation on the dividing curve (4- 35) 
of genus 3 for the points with x = 2.5 on the x-sphere. The sheets are 
identified at the fiber over -2.1404+0.44041, (-12.2492 + 2.01131, -5.1634+ 
1.35191, -4.5915 + 0.93801, -1.5405 + 0.54291). The vector of signs is s = 
(1,1,1). 




Figure 19. Solution to the DSl equation on the dividing curve (4-35) of 
genus 3 for the points a — (— 4)^-'^^ and b = (—3)^^^ att — Q. 



and we find 



A 



0.9270 -0.9270i -0.9270i 

-1. 85411 

0.92701 -0.9270 -0.92701 
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0.9270 + 0.9270i 
P^= I 

-0.9270 



0.9270 - 0.9270i 

-0.9270 + 0.9270i 

-0.9270i 





0.9270 - 0.9270i 

-0.9270i 



The following matrices A, S, C, D provide a solution of (4.6)-(4.9): 





C 




D 




lit 

1 oJ • 



which leads to the characteristic: 6 - 

To construct a solution of the DS2~ equation on the Fermat curve, we choose the points 
a = ( — 1.5 + i)^-'^) and b = (—1.5 — ip^\ The resulting solution for the choice d = can be 
seen in Fig. [20| 




Figure 20. Solution to the DS2 equation on the Fermat curve (4-36) of 
genus 3 for the points a — (—1.5 + i)^-'^) and b = (—1.5 — 1)^^^ att — 0. 



5. Conclusion 

In this paper we have presented the state of the art of the numerical evaluation of 
solutions to integrable equations in terms of multi-dimensional theta functions associated 
to real Riemann surfaces by using an approach via real algebraic curves. It was shown that 
real hyperelliptic curves parametrized by the branch points can be treated with machine 
precision for a wide range of the parameters. Even almost degenerate situations where the 
branch points coincide pairwise can be handled as long as at least one cut stays finite. This 
approach to real hyperelliptic curves [Ml [T5j is being generalized to arbitrary hyperelliptic 
curves. 

As discussed in [16] , the main difficulty for general algebraic curves is the correct numer- 
ical identification of the branch points. The case of degenerations for given branch points 
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has not yet been studied numerically, but is planned for the future. In what concerns the 



solutions (2.22) to n-NLS^ and similar solutions to the DS and the Kadomtsev-Petviashvili 
equations, the main problem in the context of real Riemann surfaces is to find the symplec- 
tic transformation leading to the homology basis introduced in [30j, for which the solutions 
of the studied equations, with regularity conditions, can be conveniently formulated. This 
problem has been reduced to find a single g x ^-matrix for given periods and real ovals, the 
latter encoded by the matrix H. For M-curves, where the matrix H vanishes, a general for- 



mula for the characteristic (4.17)-(4.19) could be given. In the general case, an algorithm 
along the lines indicated in the previous section to find the transformation will be based on 
a sufficiently general ansatz for one of the matrices entering the symplectic transformation 
which is the subject of future work. 
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